Method and system for vorticle fluid simulation

ABSTRACT

The disclosure provides an approach for animating gases. A dynamic model is employed that accounts for stretching of gas vorticles in a stable manner, handles isolated particles and buoyancy, permits deformable boundaries of objects the gas flows past, and accounts for vortex shedding. The model models stretching of vorticity by applying a vector at the center of a stretched vorticle. High frequency eddies resulting from stretching may be filtered by unstretching the vorticle while preserving mean energy and enstrophy. To model boundary pressure, a boundary may be imposed by embedding into the gas the surface boundary and setting boundary conditions based on velocity of the boundary and the Green&#39;s function of the Laplacian. For computational efficiency, a vorticle cutoff proportional to a vorticle&#39;s size may be imposed. Vorticles determined to be similar based on a predefined criteria and distance threshold may be fused.

BACKGROUND

Field of the Invention

This disclosure provides techniques for rendering images. Morespecifically, this disclosure presents techniques for simulating andrendering a gas using vorticles.

Description of the Related Art

The simulation of gases in three dimensions has been used for visualeffects, such as those in computer animations. Traditional techniquesfor simulating gases typically use voxel grids and solve pressure insuch grids. However, the animation of gases can be complex andinefficient using such traditional techniques. In addition, controllingthe motion of gas to meet the artistic needs of an animator is oftenchallenging, as gases usually do not have a well-defined surface and canvary dramatically over time.

SUMMARY

One embodiment provides a computer implemented method for renderinganimation frames depicting gaseous matter. The method generallyincludes, for each of a plurality of time steps: creating vorticles viaat least one of vortex shedding, buoyancy, and emitting vorticles;determining a harmonic field which makes flow of the gas an idealnonviscous flow; determining a velocity field which is a sum of theharmonic field and a field induced by one or more of the vorticles onboundary points placed on one or more moving objects; advecting visualparticles, density particles, and the vorticles using the velocityfield; and rendering the visual particles in an image frame.

Further embodiments include a non-transitory computer-readable storagemedium storing instructions that when executed by a computer systemcause the computer system to perform the method set forth above, and acomputer system programmed to carry out the method set forth above.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawing(s) will be provided by the Office upon request and paymentof the necessary fee.

So that the manner in which the above recited aspects are attained andcan be understood in detail, a more particular description ofembodiments of the invention, briefly summarized above, may be had byreference to the appended drawings.

It is to be noted, however, that the appended drawings illustrate onlytypical embodiments of this invention and are therefore not to beconsidered limiting of its scope, for the invention may admit to otherequally effective embodiments.

FIG. 1 illustrates an example rendering of a gas, according to anembodiment.

FIG. 2 illustrates example velocity fields induced by a vorticle,according to an embodiment.

FIG. 3 illustrates example vorticles being stretched and unstretched,according to an embodiment.

FIG. 4 illustrates an example density particle emitting a ring ofvorticles, according to an embodiment.

FIG. 5 illustrates an example harmonic field, according to anembodiment.

FIG. 6 illustrates an example of vorticle shedding, according to anembodiment.

FIG. 7 illustrates a method for rendering a gas, according to anembodiment.

FIG. 8 illustrates a system in which an embodiment may be implemented.

DETAILED DESCRIPTION

This disclosure presents techniques for efficient and controllableanimation of gases. In contrast to traditional voxel grid basedsimulation, techniques disclosed herein employ a Lagrangianvorticle-based simulation. Generally, as used herein a vorticle is avorticity particle that produces paddle-wheel like velocity, and eachvorticle carries a volume of vorticity. A novel dynamic model isintroduced that accounts for stretching of gas vorticles in a stablemanner, handles isolated particles and buoyancy, permits deformableboundaries of objects the gas flows past, and accounts for vortexshedding in which vorticles are emitted. Vorticles are not just carriedby flow, but typically also need to be squashed and stretched (elongatedalong one axis and scaled down along other axes). One embodimentsquashes and stretches vorticles in a manner that preserves the meanvelocity and vorticity of flow. In such an embodiment, the dynamic modelmay model stretching of vorticity by applying a vector at the center ofa stretched vorticle. High frequency eddies resulting from thestretching may be filtered by unstretching the stretchable vorticle in amanner that preserves both the mean energy and enstrophy of thestretchable vorticle.

To model boundary pressure, the dynamic model imposes a boundary byembedding into the gas the surface boundary and setting boundaryconditions based on a determined velocity of the boundary and theGreen's function of the Laplacian. Deformable boundaries are formulatedin a manner that only requires dot products, rather than solving linearsystems around a colliding object as in traditional techniques. In oneembodiment, the dynamic model may also impose a vorticle cutoffproportional to a vorticle's size. Doing so may improve computationalefficiency during the rendering process. Such a vorticle cutoff providesan adjustable trade-off between the fall-off's physically basedproperties and spatially localized computations. In still anotherembodiment, vorticles that are determined to be similar based on apredefined criteria and distance threshold are fused, which alsoimproves computational efficiency.

Referring now to FIG. 1, an example rendering of a gas is depicted,according to one embodiment. Panel A shows an object 101 moving througha vertical current of gas, which is simulated using a number (e.g.,2400) of vorticles 102 _(i). To render the gas, a simulation applicationemits particles having densities into the flow of the gas, the densitiesare advected, and particle densities and velocities are splatted into avolume. Advection refers to transport by a fluid due to the fluid's bulkmotion. The resulting rendering is depicted in panel B. To improvecomputational efficiency, vorticles which are sufficiently similar andclose may be fused, certain vorticles such as those outside a viewingfrustrum may be deleted, and vorticles may also be attenuated.

In one embodiment, the equation of motion for rendering the gas may beobtained by applying the curl operator to both sides of the followingform of the Navier-Stokes Equation of a viscous incompressible Newtonianfluid, dividing both sides by ρ, and replacing ∇ρ/ρ with ∇ log(ρ)

$\begin{matrix}{{\rho\frac{d\overset{\rightarrow}{v}}{dt}} = {{\mu\;{\nabla^{2}\overset{\rightarrow}{v}}} + {\rho\overset{\rightarrow}{F}} - {{\nabla p}.}}} & (1)\end{matrix}$The following equation may then be obtained to solve for motion, wherethe vorticity {right arrow over (ω)} is defined as the curl of thevelocity {right arrow over (v)}

$\begin{matrix}{\frac{d\overset{\rightarrow}{\omega}}{dt} = {{\left( {\overset{\rightarrow}{\omega} \cdot \nabla} \right)\overset{\rightarrow}{v}} + {\frac{\mu}{\rho}{\nabla^{2}\overset{\rightarrow}{\omega}}} + {{\nabla{\log(\rho)}} \times \left( {\overset{\rightarrow}{F} - \frac{d\overset{\rightarrow}{v}}{dt}} \right)} + {\nabla{\times {\overset{\rightarrow}{F}.}}}}} & (2)\end{matrix}$Equation (2) means that the vorticity {right arrow over (ω)} evolvesover time by advecting a Lagrangian frame of reference (a particle) thatcarries the vorticity {right arrow over (ω)}, as well as stretching{right arrow over (ω)} according to the velocity {right arrow over (v)},with dynamic viscosity μ, buoyancy and boundary interaction specified bydensity ρ, and external forces {right arrow over (F)} given by

$\begin{matrix}{{\overset{\rightarrow}{F}(p)} = \left\{ {\begin{matrix}{\overset{\rightarrow}{g} + \overset{\rightarrow}{e}} & {{if}\mspace{14mu} p\mspace{14mu}{outside}\mspace{14mu}{solid}\mspace{14mu}{objects}} \\\overset{\rightarrow}{f} & {otherwise}\end{matrix},} \right.} & (3)\end{matrix}$where {right arrow over (g)} is the constant for gravity, {right arrowover (f)} is the acceleration at the objects' boundaries suitable fordeformable objects, and {right arrow over (e)} represents user definedexternal forces. Density is assumed to be strictly greater than 0. Inaddition, the velocity {right arrow over (v)} that is needed foradvection may be obtained from the vorticity {right arrow over (ω)} byinverting the curl operator with the Biot-Savart law and an irrotational(i.e., not rotating) and solenoidal field {right arrow over (h)}

$\begin{matrix}{{{\overset{\rightarrow}{u}(p)} = {\frac{1}{4\pi}{\int_{x \in {\mathbb{R}}^{3}}{\frac{{\overset{\rightarrow}{\omega}(x)} \times \left( {p - x} \right)}{{{p - x}}^{3}}\ {dx}}}}}{\overset{\rightarrow}{v} = {\overset{\rightarrow}{u} + {\overset{\rightarrow}{h}.}}}} & (4)\end{matrix}$Equation (4) means that the flow {right arrow over (v)} is the sum ofthe velocity induced by a continuum of rotations of center x, axis{right arrow over (ω)} and angle ∥{right arrow over (ω)}/∥p−x∥³, with apressure field {right arrow over (h)} that models the boundarycondition.

FIG. 2 illustrates example velocity fields induced by a vorticle,according to an embodiment. In one embodiment, a vorticle basis is usedthat accounts for stretching in a stable manner. The integral insideequation (4) may be kept, while introducing a vorticle partitioning(V_(i),{right arrow over (ω)}_(i)), where {right arrow over (ω)}_(i)denotes the vorticity field induced by vorticle i

$\begin{matrix}{{\overset{\rightarrow}{u}(p)} = {\frac{1}{4\pi}{\int_{x \in V_{i}}{\frac{{{\overset{\rightarrow}{\omega}}_{i}(x)} \times \left( {p - x} \right)}{{{p - x}}^{3}}\ {{dx}.}}}}} & (5)\end{matrix}$The singularity of equation (5) at p=x could be removed by integratinganalytically the vorticity over the partition or with a regularizingconstant. To avoid this laborious integral and avoid an arbitrarypost-simulation blurring filter size, the vorticity field may be definedinstead as the sum of the curl of the vorticle's velocity field, asshown in FIG. 2, which illustrates slices of two example velocity fieldsinduced by vorticles 201-202 in panels A and B. This is in contrast todefining vorticity as interpolated from point values, and thisguarantees a divergence free vorticity and avoids instability fromvorticity compression.

A vorticle itself is a vorticity particle defined by rotation strength{right arrow over (w)}_(i), center x_(i), and falloff φ_(i)(p). Thevelocity and vorticity fields induced by a vorticle are given by{right arrow over (v)} _(i)(p)=φ_(i) {right arrow over (w)}×(p−x_(i))  (6a){right arrow over (ω)}_(i)(p)=2φ_(i) {right arrow over (w)}_(i)+∇φ_(i)×({right arrow over (w)} _(i)×(p−x _(i))).  (6b)The falloff of a stretchable vorticle

$\phi_{i} = {s_{i}{r_{i}^{- \frac{5}{2}}\left( {1 + \frac{{{\mu_{i}\left( \frac{p - x_{i}}{r_{i}} \right)}}^{2}}{2}} \right)}^{- \frac{3}{2}}}$and its gradient ∇φ_(i) are centered at x_(i), where r_(i) is thevorticle's size, s_(i) is the stretching factor, and μ_(i) is thestretching function

${\mu_{i}(q)} = {\frac{1}{{\overset{\rightarrow}{w}}_{i}^{2}}{\left( {{{S_{i}^{- \frac{4}{5}}\left( {{\overset{\rightarrow}{w}}_{i} \cdot q} \right)}{\overset{\rightarrow}{w}}_{i}} + {s_{i}^{\frac{7}{10}}{\overset{\rightarrow}{w}}_{i} \times q \times {\overset{\rightarrow}{w}}_{i}}} \right).}}$The following properties are satisfied by φ_(i). First, φ_(i) revolvesaround {right arrow over (w)}_(i), and therefore {right arrow over(v)}_(i) is divergence-free since its magnitude is constant along thestreamlines of the rotation. Second, when the stretching factor isincreased from s_(i) to s_(i)′, the vorticity {right arrow over (ω)}_(i)is stretched by a factor

$\frac{s_{i}^{\prime}}{s_{i}}$along {right arrow over (w)}_(i) and squashed by

$\sqrt{\frac{s_{i}^{\prime}}{s_{i}}}$along any direction perpendicular to {right arrow over (w)}_(i), inaccordance with the deformation induced by an incompressible flow, andsatisfying Kelvin's circulation theorem. Third, stretching isconservative since the vorticle's mean energy E_(i) is independent ofs_(i)

$\begin{matrix}{E_{i} = {{\frac{1}{2}{\int{{\overset{\rightarrow}{v}}_{i}^{2}{dx}}}} = {\sqrt{2}\pi^{2}{{{\overset{\rightarrow}{w}}_{i}}^{2}.}}}} & (7)\end{matrix}$A stretchable vorticle is thus defined by 4 variables {x_(i), {rightarrow over (w)}_(i), r_(i), s_(i)}, and, as discussed in greater detailbelow, can further be reduced to a vorticle of 3 variables{x _(i) ,{right arrow over (w)} _(i) ,r _(i)}.  (8)The velocity field is defined by summing the velocity field of multiplevorticles, as shown in FIG. 2, panel C. In panel C, three vorticles203-205 orthogonal to a plane are shown, and the velocity field is thesum of the velocity fields of the individual vorticles.

FIG. 3 illustrates example vorticles being stretched and unstretched,according to an embodiment. As discussed, in vortex simulators,vorticles are not only carried by flow, but also need to be squashed andstretched. In one embodiment, discussed in greater detail below, suchsquashing and stretching is performed in a manner that preserves themean velocity and vorticity of the flow. In particular, vorticles aremodified so that they evolve while maintaining constant enstrophy andmean energy. Illustratively, the vorticle 301 in panel A is stretchedinto the vorticle 302 in panel B. The vorticle 302 in panel B is thenresampled in an unstretched vorticle 303 shown in panel C, with the samemean energy and enstrophy as the vorticle 302 in panel B.

The term ({right arrow over (ω)}·∇){right arrow over (v)} in equation(2) models the stretching of vorticity. Stretching at a particular pointx_(i) may be measured by applying the velocity gradient to {right arrowover (ω)}_(i)(x_(i))=2r_(i) ^(−5/2){right arrow over (w)}_(i), theself-induced vorticity at the center of an unstretched vorticle. Doingso gives

$\begin{matrix}{{\frac{d{\overset{\rightarrow}{\omega}}_{i}}{dt}\left( x_{i} \right)} = {\sum\limits_{j}{{\overset{\rightarrow}{w}}_{j} \times \left( {{{{\nabla{\phi_{j}\left( x_{i} \right)}} \cdot {{\overset{\rightarrow}{\omega}}_{i}\left( x_{i} \right)}}\left( {x_{i} - x_{j}} \right)} + {{\phi_{j}\left( x_{i} \right)}{{{\overset{\rightarrow}{\omega}}_{i}\left( x_{i} \right)}.}}} \right.}}} & (9)\end{matrix}$Stretching produces a rotation of {right arrow over (w)}_(i) and scalingof s_(i).

${{Let}\mspace{14mu}{\overset{\rightarrow}{\omega}}_{i}^{\prime}} = {{\overset{\rightarrow}{\omega}}_{i} + {D_{t}{\frac{d{\overset{\rightarrow}{\omega}}_{i}}{dt}.}}}$The new rotation strength and stretch factors are then:

$\begin{matrix}{{\overset{\rightarrow}{w}}_{i}^{\prime} = \frac{{{\overset{\rightarrow}{w}}_{i}}{\overset{\rightarrow}{\omega}}_{i}^{\prime}}{{\overset{\rightarrow}{\omega}}_{i}^{\prime}}} & (10) \\{s_{i}^{\prime} = {\frac{{\overset{\rightarrow}{\omega}}_{i}^{\prime}}{{\overset{\rightarrow}{\omega}}_{i}}.}} & (11)\end{matrix}$The accumulation of stretching introduces increasingly high frequencyvelocities by transferring large scale eddies to smaller scale eddies.Although diffusion may filter eddies over long enough periods of time,instability may not be an option and high frequency eddies may need tobe filtered explicitly in a predictable manner. Such filtering may beaccomplished by unstretching the stretchable vorticle, as shown in panelC of FIG. 3, in a manner that preserves both E_(i) and the enstrophyΩ_(i) of the unstretched vorticle

$\begin{matrix}{\Omega_{i} = {{\frac{1}{2}{\int{{\overset{\rightarrow}{\omega}}_{i}^{2}{dx}}}} = {\frac{3{\pi^{2}\left( {1 + {4s_{i}^{\prime^{3}}}} \right)}}{16\sqrt{2}r_{i}^{2}s_{i}^{\prime^{8/5}}}{{{\overset{\rightarrow}{w}}_{i}}^{2}.}}}} & (12)\end{matrix}$Preserving E_(i) is trivial, as E_(i) is independent of s_(i) and r_(i)in equation (7). To preserve Ω_(i), r_(i) may be adjusted to a new sizer_(i)′

$\begin{matrix}{r_{i}^{\prime} = {r_{i}s_{i}^{\prime^{4/5}}{\sqrt{\frac{5}{1 + {4s_{i}^{\prime^{3}}}}}.}}} & (13)\end{matrix}$With equation (13), the stretching factor can be restored to 1, asillustrated in panel C of FIG. 3. It can be verified that swapping{r_(i),s_(i)′} for {r_(i)′,1} preserves Ω_(i). This step is a resamplingstep, where the same vorticle locations may be used. Note thatresampling introduces an error, especially when squashing the vorticle,i.e., when the squashing factor is below 1. If substeps are taken,unstretching may be performed on full frames to avoid overfiltering. Inone embodiment, limit resolutions may be set with a lower thresholdr_(min) and upper threshold r_(max) on the vorticle size r_(i). Such alimit resolution loses enstrophy, but does not lose energy because ofequation (7).

Equations (10)-(11) and (13) provide a way to apply stretching to avorticle. The falloff and falloff gradient for an unstretched vorticleare

$\begin{matrix}{{\phi_{i} = {\sqrt{r_{i}}\left( {r_{i}^{2} + \frac{{{p - x_{i}}}^{2}}{2}} \right)^{- \frac{3}{2}}}}{{\nabla\phi_{i}} = {{- \frac{3}{2}}\sqrt{r_{i}}\left( {r_{i}^{2} + \frac{{{p - x_{i}}}^{2}}{2}} \right)^{- \frac{5}{2}}{\left( {p - x_{i}} \right).}}}} & (14)\end{matrix}$When a vorticle becomes too small and approaches the fluid's Kolmogorovlength, viscous forces dominate, and the vorticle strength is dissipatedwith {right arrow over (w)}_(i)′=k{right arrow over (w)}_(i).

The dynamic model of equation (2) further includes a viscosity modelthat handles isolated particles. The viscosity model generally reducesthe rotational strength of vorticles by an amount proportional to ameasure of the number of surrounding vorticles or the surroundingemptiness. That is, isolated vorticles slow down by interacting withemptiness. The term

$\frac{\mu}{\rho}{\nabla^{2}\overset{\rightarrow}{\omega}}$in equation (2) represents the diffusion of vorticity and modelsviscosity. In one embodiment,

$\frac{\mu}{\rho}$may be approximated with a constant kinematic viscosity ν, and theviscous model may be derived from the modified particle strengthexchange (PSE) method, which normalizes the discrete integral to avoid ablow up. A term may also be added for handling isolated particles tomodel effectively the leakage of vorticity into the region of space withno vorticles. The PSE method is obtained by a Taylor expansion of {rightarrow over (ω)}, reduced after multiplication with a normalizedregularization function η_(ε). The result of this term is similar toartificial damping, but within the scope of the diffusion

$\begin{matrix}{{\frac{d{\overset{\rightarrow}{w}}_{i}}{dt} = {\frac{v}{\varepsilon^{2}}\left( {{\left( {1 - \alpha_{i}} \right)\frac{\sum\limits_{j}{V_{j}{\eta_{\varepsilon}\left( {x_{j} - x_{i}} \right)}\left( {{\overset{\rightarrow}{w}}_{j} - {\overset{\rightarrow}{w}}_{i}} \right)}}{\sum\limits_{j}{V_{j}{\eta_{\varepsilon}\left( {x_{j} - x_{i}} \right)}}}} - {\alpha_{i}{\overset{\rightarrow}{w}}_{i}}} \right)}}{{\alpha_{i} = {\prod\limits_{j \neq i}\left( {1 - \frac{\eta_{\varepsilon}\left( {x_{j} - x_{i}} \right)}{\eta_{\varepsilon}(0)}} \right)}},}} & (15)\end{matrix}$where α_(i) is a measure of the isolation of particle i, η_(ε) is theGaussian PSE kernel, and V_(i) is the volume associated with vorticle i

$\begin{matrix}{{{\eta_{\varepsilon}(x)} = {\frac{1}{\sqrt{2}\varepsilon^{3}\pi^{3/2}}{\exp\left( {- \frac{{x}^{2}}{2\varepsilon^{2}}} \right)}}}{V_{i} = {{\int{\phi_{i}{dx}}} = {2\sqrt{2}\pi^{2}\sqrt{r_{i}}{s_{i}^{2/5}.}}}}} & (16)\end{matrix}$By using ε=√{square root over (νD_(t))}, the viscosity is more cheaplycalculated for low ν and small time steps. This result can also beinterpreted as a convolution with the heat kernel. Contrary to howviscosity was previously handled based on fluid mechanics, in whichsimulations were performed even in areas that were not visible, theapproach discussed herein allows viscosity to leak into nothingness andvanish naturally instead of requiring a boundary that maintainsviscosity around the data even in areas that are not visible.

FIG. 4 illustrates an example density particle emitting a ring ofvorticles, according to an embodiment. To simulate buoyancy, vorticles(e.g., vorticles 401-406) are generated in rings around densityparticles carried by the flow, thereby creating a gust of wind up ordown. The density particles can be user-configured to alter thebuoyancy. When the position p is outside of solid objects, the term

${{\nabla{\log(\rho)}} \times \left( {\overset{\rightarrow}{F} - \frac{d\overset{\rightarrow}{v}}{dt}} \right)} + {\nabla{\times \overset{\rightarrow}{F}}}$in equation (2) reduces to

${\nabla\log}(\rho) \times \left( {\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{v}}{dt}} \right)$and models the vorticity induced by buoyancy. New vorticles 401-406 arethen produced from the density field ρ releasing potential energy, asshown in FIG. 4. Let ρ be defined with a set of particles carryingdensity and an ambient density ρ_(A)>0 so that the total density field ρis strictly greater than 0, as required by equation (2)ρ=ρ_(A)+Σρ_(j).  (17)Here, the subscript j is used to denote density particles, as opposed tosubscript i for vorticles. The falloff ρ_(j) of a density particle j iscentered at x_(j) and may be defined in local coordinates q_(j)=p−x_(j)as

$\begin{matrix}{\rho_{j} = {m_{j}{\frac{{\exp\left( \left( {1 + {\kappa_{1}\frac{q_{j} \cdot q_{j}}{2r_{j}^{2}}}} \right)^{- \kappa_{2}} \right)} - 1}{e - 1}.}}} & (18)\end{matrix}$where m_(j) is the multiplier of the particle density field and κ_(i)are fitting constants. As shown in FIG. 4, the newly induced vorticlesare located along a ring of diameter r_(j) perpendicular to

$\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - {\frac{d\overset{\rightarrow}{v}}{dt}.}$Instead of advecting an additional filament representation, the ring maybe discretized with n new vorticles, equidistant for simplicity, andwhere n=2 in practice. Let the orthogonal unit vectors {right arrow over(a)} and {right arrow over (b)} be defined such that the cross product{right arrow over (a)}×{right arrow over (b)} has the direction of

$\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - {\frac{d\overset{\rightarrow}{v}}{dt}.}$Given n randomly selected samples α_(j), the new vorticles are, in thedensity particle's local coordinates,

$\begin{matrix}{\mspace{79mu}{{x_{\alpha_{j}} = {x_{j} + {\frac{r_{j}}{2}\left( {{{\cos\left( \alpha_{j} \right)}\overset{\rightarrow}{a}} + {{\sin\left( \alpha_{j} \right)}\overset{\rightarrow}{b}}} \right)}}}{{\overset{\rightarrow}{w}}_{\alpha_{j}} = {r_{j}\sqrt{r_{j}}\kappa_{0}\frac{\pi}{n}{{\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{v}}{dt}}}\frac{{\rho_{j}\left( x_{\alpha} \right)} + \frac{m_{j}}{e - 1}}{\rho_{j}\left( x_{\alpha} \right)}\left( {{{{–sin}\left( \alpha_{j} \right)}\overset{\rightarrow}{a}} + {{\cos\left( \alpha_{j} \right)}\overset{\rightarrow}{b}}} \right)}}\mspace{79mu}{r_{\alpha_{j}} = {r_{j}.}}}} & (19)\end{matrix}$

FIG. 5 illustrates an example harmonic field, according to anembodiment. In one embodiment, the pressure field for boundaries is aharmonic field such as the field surrounding boundaries such as object503, the initial position of which is shown in panel A. As shown inpanel B, the harmonic field warps space from 501 to 502 in anincompressible and irrotational manner, with a slip boundary. Theharmonic field is an induced velocity field derived from the velocity atthe boundary of objects that affects vorticles around it. Conversely,new vorticles that are created may induce flows that change the harmonicfield. The harmonic field is denoted herein by {right arrow over (h)}.

In equation (4), the harmonic field {right arrow over (h)} cancels thefield induced by vorticles without adding vorticity. To define {rightarrow over (h)},

³ may be split with a surface δΩ enclosing volume Ω with normal {rightarrow over (n)} pointing outside, and a volume Ω^(c) the complementaryof Ω. It can be shown that the Neumann boundary condition then defines{right arrow over (h)} restricted to Ω^(c) as{right arrow over (h)} _(Ω) _(c) (p)=

_(δΩ) {right arrow over (n)}·({right arrow over (v)} _(Ω) _(c)(x)−{right arrow over (u)}(x))∇G(p,x)dx.  (20)Here,

${G\left( {p,x} \right)} = {- \frac{1}{4\pi{{p - x}}}}$is function of the Laplacian, and {right arrow over (v)}_(Ω) _(c) is thevelocity of the boundary. Using n panels of size a_(i) and centroidx_(i), {right arrow over (h)} may be discretized using the falloffs{right arrow over (h)} _(Ω) _(c) (p)=Σ_(i) a _(i) {right arrow over (n)}_(i)·({right arrow over (v)} _(Ω) _(c) (x _(i))−{right arrow over (u)}(x_(i)))∇G(p,x _(i)).  (21)To avoid the singularities of G when p approaches the boundary samplesx_(i), the pathlines of the monopole based on a deformer may be usedinstead of equation (21). In such a case, for a point outside of theboundary,

$\begin{matrix}{{{{\overset{->}{h}}_{\Omega^{c}}(p)} = {{\frac{1}{D_{t}}{\sum\limits_{i}{\zeta\left( {{p - x_{i}},{D_{t}a_{t}{{\overset{->}{n}}_{i} \cdot \left( {{{\overset{->}{v}}_{\Omega^{c}}\left( x_{i} \right)} - {\overset{->}{u}\left( x_{i} \right)}} \right)}}} \right)}}} - \left( {p - x_{i}} \right)}},} & (22)\end{matrix}$where ζ, defined below, satisfies ζ(ζ(p,k₀),k_(i))=ζ(p,k₀+k₁),

$\begin{matrix}{{\zeta\left( {p,k} \right)} = \left\{ {{\begin{matrix}\left( {1 + \frac{3k}{4\pi{p}^{3}}} \right)^{1/3} & {{p\mspace{14mu}{if}\mspace{14mu}{r(k)}} < {p}} \\0 & {otherwise}\end{matrix}{r(k)}} = {\left( {{\max\left( {{- k},0} \right)}{3/4}\pi} \right)^{1/3}.}} \right.} & (23)\end{matrix}$This provides a geometric insight: a boundary opposing the flow is akinto an insertion and removal of volume at the boundary proportional tothe boundary's opposition to the flow. The discretization of {rightarrow over (h)} is stable, but more accurate away from the boundary thannear the boundary. To remedy the problem of accuracy, the solution{right arrow over (h)}_(δΩ) may be defined on the boundary. Since {rightarrow over (n)}·{right arrow over (h)}_(δΩ)={right arrow over(n)}·({right arrow over (v)}_(Ω) _(c) −{right arrow over (u)}) and{right arrow over (h)}_(δΩ) is aligned with the normal{right arrow over (h)} _(δΩ)(p)=({right arrow over (n)}·({right arrowover (v)} _(Ω) _(c) (p)−{right arrow over (u)}(p))){right arrow over(n)}.  (24)In addition, if a point p enters the boundary, the point is pushed outto the nearest position on the surface{right arrow over (h)} _(Ω)(p)=argmin_(xεδΩ)(∥p−x∥).  (25)

Equations (22) and (24)-(25) are assembled to construct the fulldefinition of {right arrow over (h)}, by blending {right arrow over(h)}_(Ω) _(c) and {right arrow over (h)}_(δΩ) with a smoothstep functionbased on the distance {circumflex over (r)} between the samples on theboundary

$\begin{matrix}{{\overset{->}{h}(p)} = \left\{ {\begin{matrix}{{\overset{->}{h}}_{\Omega}(p)} & {{{if}\mspace{14mu} p} \in \Omega} \\{{mix}\left( {{{\overset{->}{h}}_{\delta\Omega}(p)},{{\overset{->}{h}}_{\Omega^{c}}(p)}} \right)} & {otherwise}\end{matrix},} \right.} & (26)\end{matrix}$where mix({right arrow over (a)},{right arrow over (b)})={right arrowover (a)}+smoothstep

${\left( \frac{d}{\hat{r}} \right)\left( {\overset{->}{b} - \overset{->}{a}} \right)},$given the distance d from p to δΩ.

FIG. 6 illustrates an example of vorticle shedding, according to anembodiment. As shown in panel A, a constant flow moves around a staticpole 610, and the surface of the pole 610 sheds vorticles. In general,an object that accelerates in a flow will stick to the flow and shedvorticles into the flow. Panel B shows that a slice of the vorticityinduced by the vorticles reveals the emergent behavior of a von Kármánvortex street.

When the position p is at the boundary of a solid object, the term

${{\nabla{\log(\rho)}} \times \left( {\overset{->}{F} - \frac{d\overset{->}{v}}{dt}} \right)} + {\nabla{\times \overset{->}{F}}}$in equation (2) measures the change of vorticity at the moving object'sboundary. This vorticity spreads into the flow by diffusion proportionalto viscosity coefficient ν introduced above. It can be shown that thesurface vorticity that satisfies the boundary condition discussed aboveis

$\left( {\overset{->}{f} - \overset{->}{e} - \frac{d\overset{->}{v}}{dt}} \right) \times {\overset{->}{n}.}$The vorticity shedding is then the solution to a differential equationwith boundary condition

$\begin{matrix}\left\{ {\begin{matrix}{\frac{d\overset{->}{\omega}}{dt} = {\upsilon{\nabla^{2}\overset{->}{\omega}}}} \\{\overset{->}{\omega} = {{{{}_{p \in \Omega}^{}{}_{}^{}}\left( {\overset{->}{f} - \overset{->}{e} - \frac{d\overset{->}{v}}{dt}} \right)} \times \overset{->}{n}}}\end{matrix}.} \right. & (27)\end{matrix}$Equation (27) can be solved by shedding vorticles. The surface may bedivided into n panels of area a_(i) and centroid x_(i), and emit perpanel a vorticle that approximates the heat kernel

$\begin{matrix}{{{\overset{->}{w}}_{i} = {a_{i}\frac{r^{5/2}}{8\pi\; D_{t}\upsilon}\left( {\overset{->}{f} - \overset{->}{e} - \frac{d\overset{->}{v}}{dt}} \right) \times \overset{->}{n}}}{r_{i} = {\sqrt{2}{\sqrt{D_{t}\upsilon}.}}}} & (28)\end{matrix}$

${\left( {\overset{->}{f} - \overset{->}{e} - \frac{d\overset{->}{v}}{dt}} \right) \times \overset{->}{n}}$can be used as a probability distribution function to create samples atthe locations that most affect the flow. Note that {right arrow over(f)} is the surface acceleration, as opposed to the velocity. To computethe acceleration, the equation (4) at the previous time may be stored onthe surface, and acceleration may be computed as

$\frac{d\overset{->}{v}}{dt} \simeq {\frac{{\overset{->}{v}(t)} - {\overset{->}{v}\left( {t - D_{t}} \right)}}{D_{t}}.}$As shown in FIG. 6, this produces the expected vorticle sheddingbehavior.

FIG. 7 illustrates a method 700 for rendering a gas, according to anembodiment. The steps of the method 700 may be repeated for a number oftime steps to simulate the gas, but only a single time step is describedwith respect to FIG. 7 for conciseness. It is assumed that, prior to thegas being rendered, a user has set up the rendering simulation bydefining, e.g., object(s) in a scene, colliders, gravity, gust(s) ofwind, and emitters of density, among other things. In one embodiment,three types of particles may be defined by the user: vorticles forrotations, density particles for buoyancy, and visual particles that areactually rendered to an image and controlled by the vorticles, densityparticles, and object boundaries during simulation of the gas. Thevorticles may each have an axis around which rotation is happening, astrength of the rotation, and a size representing the footprint of thevorticle in space; the density particles may have mass and work inconnection with gravity; and the visual particles may have density. Insuch a case, a user may place vorticles if the user wishes to create awind pushing one way or another, the user may place density particles ina scene if the user wishes the gas to float or sink, and the user mayplace visual particles to be able to see the simulation result. Inaddition, the user may define deformable boundaries of objects in thescene which are sampled as colliders during the simulation.

In another embodiment, the number of samples of the colliders, shedding,sources, buoyancy, and density may be controlled independently, andusers are able control a flow by modifying existing vorticles withexternal forces {right arrow over (e)} or via the harmonic field {rightarrow over (h)}, or by creating new vorticles. For example, a turbulentfield may be created by scattering vorticles with random parameters, agust of wind may be created by placing vorticles aligned with thetangent of a ring perpendicular to the direction of the desired wind, aninvisible collider may be moved to warp space, and the amount ofshedding of vorticles and buoyancy can be artificially dialed up ordown. Additional dials that may be tunable per vorticle in someembodiments include overshoot, spin, and artificial damping. Overshootand undershoot refer to, for advection of visible particles, the scalinghigher or lower of velocity to create drag or extra swirliness. Althoughphysically implausible, this can help create styles. Spin refers to theaxis of vorticles being aligned with vector fields to modify the fluidmotion globally. Artificial damping refers to strength of vorticlesbeing artificially reduced using, e.g., a damping coefficient that isstored per vorticle.

At step 710, the simulation application scatters point samples onboundaries of objects using the probability distribution discussed abovewith respect to FIG. 6. As discussed, the surface of objects are dividedinto areas that emit per area a monopole that approximates a surfacepixel, and the probability distribution

${\left( {\overset{->}{f} - \overset{->}{e} - \frac{d\overset{->}{v}}{dt}} \right) \times \overset{->}{n}}$may be used to create samples at locations that most affect the flow inone embodiment.

At step 720, the simulation application creates vorticles via one or acombination of vortex shedding, buoyancy, and emitting vorticles definedby a user. In one embodiment, vortex shedding may be defined by equation(28), discussed above. For example, an object accelerating through aflow may shed vorticles. New vorticles may also be created frombuoyancy, as defined by equation (19), discussed above. In addition, newvorticles may be emitted as defined by the user. For example, the usermay scatter a uniform distribution of random vorticles to produce aninitial condition, and vorticles may be emitted on a circle withtangential strength, producing a source of wind.

At step 730, the simulation application computes a harmonic field whichmakes flow of the gas an ideal nonviscous flow. In one embodiment, theharmonic field {right arrow over (h)} induced by colliders on visualparticles, density particles, and vorticles may be computed usingequation (26), discussed above.

At step 740, the simulation application computes a velocity field thatis the sum of the harmonic field and the field induced by vorticles onboundary points placed on one or more of the object(s) in the scene. Inone embodiment, the velocity {right arrow over (u)} induced by vorticleson visual particles that are actually seen by the user, densityparticles, and vorticles may be computed using equation (4), discussedabove.

Then, at step 750, the simulation application advects visual particles,density particles, and vorticles of the gas using the velocity field. Ingeneral, the advecting process stretches or squashes (a squash being areverse stretch) vorticles, thereby modifying the shape of thevorticles' footprints as well as potentially changing the position andaxis of rotation of the vorticles. As discussed, stretched vorticles mayfurther be unstretched by measuring a new footprint for the vorticlesuch that a new vorticle that is round has the same enstrophy and meanenergy in a next image frame.

In one embodiment, the simulation application may apply the displacementD_(t)({right arrow over (u)}+{right arrow over (h)}) to visualparticles, density particles, and vorticles during the advectingprocess. Several integration schemes may be used to advect particlesalong an induced velocity field. A simple method which is stable forlarge steps relies on circular path lines of individual vorticles andreplaces equation (6a), above, with

$\begin{matrix}{{{\overset{->}{k}}_{w} = {D_{t}{\phi_{i}\left( {p - x_{i}} \right)}{\overset{->}{w}}_{i}}}{{\overset{->}{k}}_{x} = {{\overset{->}{k}}_{w} \times \left( {p - x_{i}} \right)}}{{{\overset{->}{v}}_{i}(p)} = {\frac{1}{D_{t}}{\left( {{\frac{1 - {\cos\left( {{\overset{->}{k}}_{w}} \right)}}{{\overset{->}{k}}_{w}^{2}}{\overset{->}{k}}_{w} \times {\overset{->}{k}}_{x}} + {\frac{\sin\left( {{\overset{->}{k}}_{w}} \right)}{{\overset{->}{k}}_{w}}{\overset{->}{k}}_{x}}} \right).}}}} & (29)\end{matrix}$

At step 760, the simulation application fuses vorticles, as appropriate.Fusing vorticles that are sufficiently similar to each other andsufficiently close in distance may improve computational efficiency. Thefusing may be performed in a hierarchical manner, as is known to one ofskill in the art. In one embodiment, vorticles may be considered similarenough when ∥p_(i)−p_(j)∥ and |r_(i)−r_(j)| are smaller than a distancethreshold ε which is user-controllable. Also, the velocity induced by agroup of far away vorticles can be obtained by fusing the vorticles. Thefused vorticle is given by the following formula

$\begin{matrix}{{x = \frac{\sum\; x_{i}}{n}}{r = \left( \frac{\sum\limits_{i}r_{i}^{- \frac{5}{2}}}{\sum\limits_{i}\sqrt{r_{i}}} \right)^{{- 1}/3}}{{\overset{->}{w} = \frac{\sum\limits_{i}{\sqrt{r_{i}}{\overset{->}{w}}_{i}}}{\sqrt{r}}},}} & (30)\end{matrix}$where n is the number of vorticles in the cell and i is the vorticleindex. The above is asymptotic to the sum of vorticles and equal to thesum when the vorticles have the same radius.

Then at step 770, the simulation application attenuates and deletesvorticles, as appropriate. This step includes viscosity. Similar tofusing vorticles, attenuating and deleting vorticles may improvecomputational efficiency. The attenuation of vorticles may be based onthe diffusion of vorticles leaking vorticity to regions which contain novorticity, and the rate of attenuation may be user-controllable.Attenuated vorticles may eventually cease to spin and be deleted. Thevorticles that are deleted may include vorticles far from visualparticles, vorticles outside of a view frustrum, vorticles having lowstrength, or some combination of these.

In one embodiment, a vorticle cutoff proportional to vorticle size r_(i)may be used. Doing so reduces the algorithmic complexity of the method700 to O(N) when neighbor cell lists are used. This accelerates thenearest vorticle search while preserving an incompressible flow sincethe cutoff is constant along the streamlines of rotation. Sincevorticles are band limited, the vorticles may further be split intogroups of similar radii, and their displacements evaluated in sparsegrids with cell size proportional to the cutoff. In a particularembodiment, the vorticle cutoff may be introduced on both the falloffand its gradient in a manner that preserves smoothness. Let thevorticle's cutoff distance be defined as mr_(i). Using the followingvalues

$\begin{matrix}{{k_{0} = {4\left( {1 + {2m^{2}}} \right)r_{i}}}{k_{1} = {2 + m^{2}}}{k_{2} = {k_{1}^{2}\sqrt{k_{1}}}}{a_{0} = \frac{r_{i}}{r_{i} - {\sqrt{2}{\kappa_{0}/\kappa_{2}}}}}{a_{1} = {\sqrt{2}\frac{\kappa_{0} - {6m{{p - x_{i}}}}}{\kappa_{2}\sqrt{r_{i}}r_{i}^{3}}}}{{a_{3} = \frac{6\sqrt{2}}{\kappa_{2}\sqrt{r_{i}}r_{i}^{4}}},}} & (31)\end{matrix}$the remapped falloff and gradient are{tilde over (φ)}_(i) =a ₀(φ_(i) −a ₁)∇{tilde over (φ)}_(i)=∇φ_(i) +a ₃(p−x _(i)).  (32)This vorticle cutoff provides an adjustable tradeoff between thefalloff's physically based properties and spatially localizedcomputations, with defaults set to, e.g., 6r_(i) and r_(i)/2.

In another embodiment, vorticles may be assigned an artificial decayrate, triggered by an event controlled by a varying expression. As aresult, the vorticles may have limited lifespan. In a furtherembodiment, band-limiting may be employed in which, as stretchingdeforms vorticles outside of the range (r_(min), r_(max)), thevorticle's strength {right arrow over (w)}_(j) is reduced instead ofscaling the vorticle's radius r_(i). For shrinking small vorticles, thismodels the fluid's viscous behavior at small scales, and for expandinglarge vorticles this reduces their contribution to velocityartificially. In an additional embodiment, vorticle radii may be paged.Splitting the vorticles in groups of similar radii has been seen to leadto more efficient use of acceleration structure for distance queries,and a logarithmic scale may be used to split vorticles in such groups.In yet another embodiment, velocities may be cached. When the pointdensity is particularly high within a region of space, velocity may beevaluated at the vertices of a lattice and interpolated in-between.

At step 780, the simulation application renders the visual particles inan image frame, which is the output viewable by a user.

FIG. 8 depicts a block diagram of a system 800 in which an aspect ofthis disclosure may be implemented. As shown, the system 800 includes,without limitation, a central processing unit (CPU) 810, a networkinterface 830, an interconnect 815, a memory 860 and storage 820. Thesystem 800 may also include an I/O device interface 840 connecting I/Odevices 850 (e.g., keyboard, display and mouse devices) to the system800.

The CPU 810 retrieves and executes programming instructions stored inthe memory 860. Similarly, the CPU 810 stores and retrieves applicationdata residing in the memory 860. The interconnect 815 facilitatestransmission, such as of programming instructions and application data,between the CPU 810, I/O device interface 840, storage 820, networkinterface 830, and memory 860. CPU 810 is included to be representativeof a single CPU, multiple CPUs, a single CPU having multiple processingcores, one or more graphics processor units (GPUs), and the like, orsome combination of these. And the memory 860 is generally included tobe representative of a random access memory. The storage 820 may be adisk drive storage device. Although shown as a single unit, the storage820 may be a combination of fixed and/or removable storage devices, suchas fixed disc drives, removable memory cards or optical storage, networkattached storage (NAS), or a storage area-network (SAN). Further, system800 is included to be representative of a physical computing system aswell as virtual machine instances hosted on a set of underlying physicalcomputing systems. Further still, although shown as a single computingsystem, one of ordinary skill in the art will recognized that thecomponents of the system 800 shown in FIG. 8 may be distributed acrossmultiple computing systems connected by a data communications network.

As shown, the memory 860 includes an operating system 861 and asimulation application 862. For example, the operating system may beMicrosoft Windows®. The simulation application 862 is configured tosimulate a gas. In one embodiment, the simulation application 862 may,for a number of time steps: scatter point samples on boundaries ofobjects; create vorticles via one or a combination of vortex shedding,buoyancy, and emitting vorticles defined by a user; compute a harmonicfield which makes flow of the gas an ideal nonviscous flow; compute avelocity field that is the sum of the harmonic field and the fieldinduced by vorticles on boundary points placed on one or more of theobject(s) in the scene; advect visual particles, density particles, andthe vorticles using the velocity field; fuse vorticles, attenuatevorticles, and delete vorticles, as appropriate; and render the visualparticles in an image, as discussed above with respect to FIG. 7.Rendered images may then be stored in the storage 820.

Advantageously, techniques disclosed herein permit gases to be simulatedusing vorticles, providing a full set of features expected from a gassimulation while using only vorticles and monopoles. In contrast totraditional voxel grid based simulation, techniques disclosed hereinemploy a Lagrangian vorticle-based simulation which uses point integralsand is more scalable than traditional techniques that required computingmatrix inverses. The vorticles are defined in relation to theNavier-Stokes equation. The simulation of gases has no divergence byconstruction and avoids solving pressure completely. The simulation isalso gridless and can span domains of arbitrary size and resolution. Thesampling resolutions of density, dynamics, boundary condition andsources are decoupled from each other, which permits many options foradaptive sampling strategies. In addition, techniques disclosed hereinadvect vorticles that carry a volume of vorticity, rather than filteringa point vorticity. Doing so leads to stable stretching. Furtherstability is achieved by integrating streamlines instead of velocity.The model for boundary pressure disclosed herein does not requiresolving a linear system, and the model disclosed herein for buoyancywhich plays a role in gas motion is capable of producing lively smokesand pyroclastic effects.

Reference is made herein to embodiments of the invention. However, itshould be understood that the invention is not limited to specificdescribed embodiments. Instead, any combination of the describedfeatures and elements, whether related to different embodiments or not,is contemplated to implement and practice the invention. Furthermore,although embodiments of the invention may achieve advantages over otherpossible solutions and/or over the prior art, whether or not aparticular advantage is achieved by a given embodiment is not limitingof the invention. Thus, the aspects, features, embodiments andadvantages are merely illustrative and are not considered elements orlimitations of the appended claims except where explicitly recited in aclaim(s). Likewise, reference to “the invention” shall not be construedas a generalization of any inventive subject matter disclosed herein andshall not be considered to be an element or limitation of the appendedclaims except where explicitly recited in a claim(s).

As will be appreciated by one skilled in the art, aspects of the presentinvention may be embodied as a system, method or computer programproduct. Accordingly, aspects of the present invention may take the formof an entirely hardware embodiment, an entirely software embodiment(including firmware, resident software, micro-code, etc.) or anembodiment combining software and hardware aspects that may allgenerally be referred to herein as a “circuit,” “module” or “system.”Furthermore, aspects of the present invention may take the form of acomputer program product embodied in one or more computer readablemedium(s) having computer readable program code embodied thereon.

Any combination of one or more computer readable medium(s) may beutilized. The computer readable medium may be a computer readable signalmedium or a computer readable storage medium. A computer readablestorage medium may be, for example, but not limited to, an electronic,magnetic, optical, electromagnetic, infrared, or semiconductor system,apparatus, or device, or any suitable combination of the foregoing. Morespecific examples (a non-exhaustive list) of the computer readablestorage medium would include the following: an electrical connectionhaving one or more wires, a portable computer diskette, a hard disk, arandom access memory (RAM), a read-only memory (ROM), an erasableprogrammable read-only memory (EPROM or Flash memory), an optical fiber,a portable compact disc read-only memory (CD-ROM), an optical storagedevice, a magnetic storage device, or any suitable combination of theforegoing. In the context of this document, a computer readable storagemedium may be any tangible medium that can contain, or store a programfor use by or in connection with an instruction execution system,apparatus, or device.

A computer readable signal medium may include a propagated data signalwith computer readable program code embodied therein, for example, inbaseband or as part of a carrier wave. Such a propagated signal may takeany of a variety of forms, including, but not limited to,electro-magnetic, optical, or any suitable combination thereof. Acomputer readable signal medium may be any computer readable medium thatis not a computer readable storage medium and that can communicate,propagate, or transport a program for use by or in connection with aninstruction execution system, apparatus, or device.

Program code embodied on a computer readable medium may be transmittedusing any appropriate medium, including but not limited to wireless,wireline, optical fiber cable, RF, etc., or any suitable combination ofthe foregoing.

Computer program code for carrying out operations for aspects of thepresent invention may be written in any combination of one or moreprogramming languages, including an object oriented programming languagesuch as Java, Smalltalk, C++ or the like and conventional proceduralprogramming languages, such as the “C” programming language or similarprogramming languages. The program code may execute entirely on theuser's computer, partly on the user's computer, as a stand-alonesoftware package, partly on the user's computer and partly on a remotecomputer or entirely on the remote computer or server. In the latterscenario, the remote computer may be connected to the user's computerthrough any type of network, including a local area network (LAN) or awide area network (WAN), or the connection may be made to an externalcomputer (for example, through the Internet using an Internet ServiceProvider).

Aspects of the present invention are described herein with reference toflowchart illustrations and/or block diagrams of methods, apparatus(systems) and computer program products according to embodiments of theinvention. It will be understood that each block of the flowchartillustrations and/or block diagrams, and combinations of blocks in theflowchart illustrations and/or block diagrams, can be implemented bycomputer program instructions. These computer program instructions maybe provided to a processor of a general purpose computer, specialpurpose computer, or other programmable data processing apparatus toproduce a machine, such that the instructions, which execute via theprocessor of the computer or other programmable data processingapparatus, create means for implementing the functions/acts specified inthe flowchart and/or block diagram block or blocks.

These computer program instructions may also be stored in a computerreadable medium that can direct a computer, other programmable dataprocessing apparatus, or other devices to function in a particularmanner, such that the instructions stored in the computer readablemedium produce an article of manufacture including instructions whichimplement the function/act specified in the flowchart and/or blockdiagram block or blocks.

The computer program instructions may also be loaded onto a computer,other programmable data processing apparatus, or other devices to causea series of operational steps to be performed on the computer, otherprogrammable apparatus or other devices to produce a computerimplemented process such that the instructions which execute on thecomputer or other programmable apparatus provide processes forimplementing the functions/acts specified in the flowchart and/or blockdiagram block or blocks

While the foregoing is directed to embodiments of the present invention,other and further embodiments of the invention may be devised withoutdeparting from the basic scope thereof, and the scope thereof isdetermined by the claims that follow.

What is claimed is:
 1. A method for rendering animation frames depictinggaseous matter, comprising, for each of a plurality of time steps:creating vorticles via at least one of vortex shedding, buoyancy, andemitting vorticles; determining a harmonic field which makes flow of thegaseous matter an ideal nonviscous flow; determining a velocity fieldwhich is a sum of the harmonic field and a field induced by one or moreof the vorticles on boundary points placed on one or more movingobjects; advecting visual particles, density particles, and thevorticles using the velocity field; and rendering the visual particlesin an image frame.
 2. The method of claim 1, wherein a basis of thevorticles permits stable stretching of the vorticles.
 3. The method ofclaim 2, wherein stretching and squashing of the vorticles is performedin a manner that maintains constant enstrophy and mean energy.
 4. Themethod of claim 1, wherein the velocity field is determined using adynamic model that permits user control of at least one of externalforces, buoyancy, rigid and deformable boundaries, and viscosity.
 5. Themethod of claim 4, wherein sampling resolutions for density, dynamics,boundary conditions, and sources in the dynamic model are decoupled fromeach other.
 6. The method of claim 4, wherein the dynamic model includesa boundary pressure model in which a boundary of the gaseous matter isimposed by embedding into the gaseous matter a surface boundary andsetting boundary conditions based on velocity of the surface boundaryand a Green's function of the Laplacian.
 7. The method of claim 4,wherein the dynamic model includes a viscosity term that reducesstrength of the vorticles based on a measure of isolation of thevorticles.
 8. The method of claim 4, wherein the dynamic model includesa buoyancy term that models density particles in the gaseous matteremitting rings of vorticles with strengths corresponding to masses ofthe density particles.
 9. The method of claim 1, wherein the vortexshedding includes emitting, for each of one or more of a plurality ofpanels into which surfaces of the objects are divided, a respectivevorticle that approximates a heat kernel, based on a probabilitydistribution function.
 10. The method of claim 1, further comprising,for one or more of the time steps, at least one of: fusing two or moreof the vorticles based on distance and similarity criteria; deleting oneor more of the vorticles which are far from the visual particles,outside a view frustrum, or whose strength is less than a thresholdvalue; and attenuating one or more of the vorticles.
 11. Anon-transitory computer-readable storage medium storing a program,which, when executed by a processor performs operations for renderinganimation frames depicting gaseous matter, the operations comprising,for each of a plurality of time steps: creating vorticles via at leastone of vortex shedding, buoyancy, and emitting vorticles; determining aharmonic field which makes flow of the gaseous matter an idealnonviscous flow; determining a velocity field which is a sum of theharmonic field and a field induced by one or more of the vorticles onboundary points placed on one or more moving objects; advecting visualparticles, density particles, and the vorticles using the velocityfield; and rendering the visual particles in an image frame.
 12. Thecomputer-readable storage medium of claim 11, wherein a basis of thevorticles permits stable stretching of the vorticles.
 13. Thecomputer-readable storage medium of claim 12, wherein stretching andsquashing of the vorticles is performed in a manner that maintainsconstant enstrophy and mean energy.
 14. The computer-readable storagemedium of claim 11, wherein the velocity field is determined using adynamic model that permits user control of at least one of externalforces, buoyancy, rigid and deformable boundaries, and viscosity. 15.The computer-readable storage medium of claim 14, wherein samplingresolutions for density, dynamics, boundary conditions, and sources inthe dynamic model are decoupled from each other.
 16. Thecomputer-readable storage medium of claim 14, wherein the dynamic modelincludes a boundary pressure model in which a boundary of the gaseousmatter is imposed by embedding into the gaseous matter a surfaceboundary and setting boundary conditions based on velocity of thesurface boundary and a Green's function of the Laplacian.
 17. Thecomputer-readable storage medium of claim 14, wherein the dynamic modelincludes a viscosity term that reduces strength of the vorticles basedon a measure of isolation of the vorticles.
 18. The computer-readablestorage medium of claim 14, wherein the dynamic model includes abuoyancy term that models density particles in the gaseous matteremitting rings of vorticles with strengths corresponding to masses ofthe density particles.
 19. The computer-readable storage medium of claim11, wherein the vortex shedding includes emitting, for each of one ormore of a plurality of panels into which surfaces of the objects aredivided, a respective vorticle that approximates a heat kernel, based ona probability distribution function.
 20. The computer-readable storagemedium of claim 11, further comprising, for one or more of the timesteps, at least one of: fusing two or more of the vorticles based ondistance and similarity criteria; deleting one or more of the vorticleswhich are far from the visual particles, outside a view frustrum, orwhose strength is less than a threshold value; and attenuating one ormore of the vorticles.
 21. A system, comprising: a processor; and amemory, wherein the memory includes an application program configured toperform operations for rendering animation frames depicting gaseousmatter, the operations comprising, for each of a plurality of timesteps: creating vorticles via at least one of vortex shedding, buoyancy,and emitting vorticles, determining a harmonic field which makes flow ofthe gaseous matter an ideal nonviscous flow, determining a velocityfield which is a sum of the harmonic field and a field induced by one ormore of the vorticles on boundary points placed on one or more movingobjects, advecting visual particles, density particles, and thevorticles using the velocity field, and rendering the visual particlesin an image frame.